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ABSTRACT 

We investigate the effect of dimensionality on the transition to explosion in neutrino- 
driven core-collapse supernovae. Using parameterized hydrodynamic simulations of 
the stalled supernova shock in one-, two- (2D), and three spatial dimensions (3D), 
we systematically probe the extent to which hydrodynamic instabilities alone can tip 
the balance in favor of explosion. In particular, we focus on systems that are well 
into the regimes where the Standing Accretion Shock Instability (SASI) or neutrino- 
driven convection dominate the dynamics, and characterize the difference between 
them. We find that SASI-dominated models can explode with up to ~ 20% lower 
neutrino luminosity in 3D than in 2D, with the magnitude of this difference decreasing 
with increasing resolution. This improvement in explosion conditions is related to 
the ability of spiral modes to generate more non-radial kinetic energy than a single 
sloshing mode, increasing the size of the average shock radius, and hence generating 
better conditions for the formation of large-scale, high-entropy bubbles. In contrast, 
convection-dominated explosions show a smaller difference in their critical heating rate 
between 2D and 3D (< 8%), in agreement with previous studies. The ability of our 
numerical implementation to maintain arbitrary symmetries is quantified with a set 
of SASI-based tests. We discuss implications for the diversity of explosion paths in a 
realistic supernova environment. 

Key words: hydrodynamics - instabilities - neutrinos - nuclear reactions, nucle¬ 
osynthesis, abundances - shock waves - supernovae: general 


1 INTRODUCTION 

While the explosion mechanism of massive stars is still an 
unsolved problem, important progress has been made in t he 
last few decades on the theory side (see, e.g.. iFoglizzo et ahl 
120151 for a recent review). At present, the leading explanation 
for the majority of core-collapse supernovae is the neutrino 
mechanism, in which energy deposition by a small fraction 
of the outgoing neutrinos po wers the revival of the stalled 
shock IIBethe fe Wilsonlll985l ). State-of-the-art models have 
shown that this mechanism fails in spherical symmetry, ex¬ 
cept for a small fraction of p rogen itors at the low mass end 
for core-collapse (e.g.. |jankall2012l and referneces therein). 


Significant effort has been invested in the develop¬ 
ment of multi-dimensio nal hydrodynamic models, first in 


axisymmetry (2D; e.e 

. Hcran 

t et al. 1992!; Burrows et al. 

1995; 

Janka & Muller 

1993; 

Ott et al. 20081; Suwa et al. 

201C; 

Bruenn et al. 20131) and recently in three spatial 


dimensions (3D), with varying degrees of approxima¬ 
tion for ot h er p hysics such as neutrino transport (e.g., 
lHanke et al.l I2OI3I and references therein). These stud¬ 


ies and others have demonstrated the importance of hy¬ 
drodynamic instabilities for providing favorable conditions 
for explosion. Observationally, an asym metric explosion is 
expec ted from spectropolarimetry (e.g., IWang fe Wheelen 
12008T ). the d istribu tion of elements in supernova remnants 


of pulsars (e.g., Hobbs et al 


(e.g. JGrefenstett e et al.||2014 ), and the large proper motions 


200511 . 


Two types of instability can break the symmetry of a 
stalled shock. First, neutrino energy deposition drives buoy¬ 
ant c onvection in a layer just inside the shock (e.g., iBethel 
Il990l l. This is a local instability that generates kinetic en¬ 
ergy on relatively small spatial scales (spherical harmonic 
indices £ > 5). Second, the shock is unstable to a global 
oscillatory instability driven by a cycle of advected and 
acoustic perturbations trapped in between the shock and 
the neutri nosphere: the Standing Accretion Shock Instabil¬ 
ity (SASI; Blondin et alj|2003l; Blondin fc Mezzacap~pall2006l : 


IFoglizzo et al.ll2007l : iGuilet fc Foglizzoll2012l l. This instabil- 
ity favors larger spatial scales than convection (£ ~ 1 — 2). 
Non-axisymmetric (spiral) modes of the SASI can also lead 
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2 Rodrigo Fernandez 


to angular momentum redistribution, as found in numerical 
simu l ations llBlondin &; Mezzacappal 2007: Blondin fc Shawl 


120071 : IFernandez 2010^_ Guilei i fc Fernandez! 20141 ') and" in ex- 

periments (IFoglizzo et al. I l2012i) . 

The first generation of 2D supernova models with ad¬ 
vanced neutrino transport found strong sloshing activity in 
the shock preceding explosion, and hence the SASI was 
generally cred i ted f o r contributing to thi s success (e.g., 
IBurrows et al.l 120071 : iMarek fc Jankal 120091 ). But because 
the inverse cascade in 2D turbulence transfers kinetic en¬ 
ergy to large spatial scales, it was not clear whether these 
sloshings were in fact driven by the SAS I or just a conse- 
guence of convective activity fe.g.. IFernandez fc Thompsonl 
I2009al ). Thereafter, the first batch of 3D hydrodynamic 
core-collapse models with approximate neutrino physics 


(Fryer & Young 

2007; 

Nordhaus et al. 20ld: Hanke et al. 

2012 

: iTakiwaki et al 

20121: Couch :2013l: Dolence et al. 


zuidll . snowing instead all tne signs 01 Buoyant convection 
drivi ng the dynamics (IBurrows et al.l 120121 ; iMurphv et al.l 
120131) . These early 3D studies however focused on a small 
set of progenitors. 

The full-physics 2D study of iMiiller et al.l (120121 ) col¬ 
lapsed two progenitors (8.1 and 27M©) with very differ¬ 
ent inner density profiles (compactness), and found two dis¬ 
tinctive explosion paths: one dominated by convection and 
another dominated by the SASI at low and high accretion 
rates, respectively. A followup study of the 27Mq progenitor 
in 3D using sophisticated neutrino transport fo und episodic 
SASI activity, but no successful explosion (lHanke et al.l 
120131 ). Other studies using more approximate tra nsport did 
not f i nd a dominant SASI in exploding system s (|Ott e t al.l 
120131 : ICouch &; O’Comioi1l20l4 lAbdikamalov et all 20141) . 


The presence of SASI activity lea ds to specific pre¬ 
dictions regardin g the n eutrino sig nal (Marck ct al. 120091 : 
iLund et al .11 201 til . 120121 : iTamborra et al .1 120131. 1 20141) and 
gravitational wave emission (e.g., Murphy et al. 120091 : 
iKotake et al.l 1201 ll : IMiiller et al.l 120131) . The differences in¬ 
troduced by dimensionality in bona fide S ASI-dominated ex¬ 
plosions have not been systematically explored vet, with the 
closes t work being the par a metric mode ls of lHanke et al.l 
d2013l) and llwakami et al.l d2008l . 120141) . in which sys¬ 
tems straddle both instability regimes (convection param¬ 
eter close to the critical value) at the time when the 
shock stalls or in the initial condition. In contrast, multi¬ 
ple 3D studies have addressed convection-dominated sys¬ 
tems, finding either relatively minor differences in the 
susceptibility to explosion between 2D and 3D (e.g ., 
lHanke et al. 120121 : iDolence et al.l 120131 : lHandv et al.l 201^1 
or m o re detrimental conditions in 3D (e.g.. lHanke et aD 
20j^_ [Couch_fc O’ Connorl I20l4 iMezzaca ppa et a hl 120151 : 


Lentz et al.l 20151 : although recent work by Miillerl 20151 sug 


gests that 3D can be more favorable than 2D after the onset 
of explosion). A better understanding of SASI-dominated 
explosions is important within the search for a robust ex¬ 
plosion mechanism, because the susceptibility to the SASI 
depends not only on the compactness of the progenitor 
via the accretion rate, but also on the speed a t which the 
protoneutron star contracts dScheck et al.ll2008l) . Therefore, 
the current generation of models could be underestimating 
the prevalence of the SASI given the treatment of neutrino 


physics and the uncertainties in the stellar models and equa¬ 
tion of state (EOS). 

Here we investigate the transition to explosion using pa¬ 
rameterized hydrodynamic models that, while approximate, 
are free from uncertainties about the progenitor and the 
EOS, allowing a clean test of the hydrodynamic effects in 
2D and 3D. Model parameters can be controlled so that 
at explosion, they lie well within the regime of dominance 
of either instabi li ty. O ur approach builds on the study of 
IFernandez et alj (|2014| . hereafter Paper I), who character¬ 
ized the properties of SASI- and convection-dominated ex¬ 
plosions in 2D. The main conclusion of that study is that 
both types of explosion succeed by generating large-scale, 
high-entropy bubbles in different ways: either directly via 
large-scale shock expansions (SASI) or through the dynam¬ 
ics of convective bubble growth. Our focus here is to probe 
the differences introduced by adding the 3rd spatial dimen¬ 
sion to SASI-dominated systems, and to compare the out¬ 
come to convection-dominated explosions. 

The structure of the paper is the following. Section [5] 
describes the physical model and numerical method used, 
and the list of models evolved. Section [3] presents our re¬ 
sults, with an overview of the transition to explosion, the 
role of dimensionality and resolution, and the properties of 
non-exploding SASI-dominated systems. Section [3] summa¬ 
rizes our results and discusses the implications for our un¬ 
derstanding of core-collapse supernovae. The Appendix de¬ 
scribes the extension of the code to 3D and a set of SASI- 
based tests to quantify its accuracy. 


2 METHODS 

2.1 Physical Model and Initial Condition 

In order to study the hydrodynamic instabilities that arise 
during the stalled shock phase of core-collapse supernovae, 
we employ a non-rotating spherical accretion flow with a 
standing shock as the baseline state. This solution assumes 
constant mass accretion rate, constant neutron star radius, 
and employs a gamma-law equation of state with param¬ 
eterized charged-current neutrino source terms, including 
‘light b u lb’ irr a diatio n. This type of model was first used by 
[Blondin et al.1 (120031 ) without heating to isolate the SASI, 
and has been subsequently employed in a n umber of para- 


metric studies of the CCSNe mechanism (Foelizzo et al. 

200f 

. 2007|: Fernandez & Thompson! 2009b. d: Guilet et al. 

201Q 

: Guilet & Foelizzo) 201^). 


This parametric approach has two advantages: (1) the 
results are independent of the uncertainties in stellar mod¬ 
els, the dense matter EOS, or an incomplete treatment 
of neutrino transport; and (2) the lin ear stability proper¬ 
ties of the system are well unde rstood (IFoglizzo et al.ll2007l : 
IFernandez fc Thompsonl l2009bl ). enabling the construction 
of model sequences in well-defined regions of parameter 
space. The resulting flow in the gain region compares fa¬ 
vorably with models that include more physics (Paper I). 

The initial condition consists of a steady-state spheri¬ 
cal accretion flow incident on a neutron star of mass M nB 
and radius r,, with a standing shock at a radius r 3 . The 
flow outside of the shock is adiabatic and supersonic, with 
Mach number Mo at a fiducial radius r 0 . Inside the shock 


© 2012 RAS, MNRAS 000.HHT81 
























































































































































SASI- and Convection-Dominated CCSNe in 3D 3 


the flow is subsonic, decelerating onto the neutron star sur¬ 
face through the emission of neutrinos and geometric con¬ 
vergence. The specific rate of neutrino heating and cooling 
is parameterized as 


Q v 


B_ 

j.2 



e (s/3n,in)2 Q(Mcut-M), 


( 1 ) 


where p is the pressure and s is the entropy. The mini¬ 
mum entropy s m i n prevents runaway cooling, while the cutoff 
Mach number M cu t suppresses source terms in the upstream 
flow. 

To generate a model sequence, we first fix the cooling 
normalization in the absence of heating (B = 0). Given the 
ratio of shock to stellar radius r s /r», with r s = ro, the adi¬ 
abatic index 7 , and the constant rate of nuclear dissocia¬ 
tion at the shock g, the cooling normalization A is found by 
requiring that the radial velocity vanishes at r = r*. The 
sequence is then formed by varying the heating constant B, 
which is related to the liglitbulb neutrino luminosity by 


B ~ 0.007 L„ e ,52 7,4 V2 ’ ‘ 2 > 


where Z /„ e ,52 = L„ e /( 10 52 erg s _1 ) is the electron neutrino 
luminosity and 2 %4 = T„/(4 MeV) is the neutrinospheric 
temperature. Increasing B increases the initial shock radius 
r a from its value without heating, ro. The minimum entropy 
Smin is obtained separately for each model, by computing 
the entropy at r = r* without using the exponential cutoff 
in equation ©■ The numerical value is obtained using the 
ideal gas entro py normalized to its post-shock value (e.g., 
iFoglizzo et al.lf2007h 


s = 


1 


7 - 1 


In 


JL 

P2 



(3) 


where p 2 and P 2 are the initial post-shock density and pres¬ 
sure. The initial condition for the time-dependent simulation 
is then obtained by re-computing the solution including the 
entropy cutoff s m i n in equation JT]). 

To connect with previous work, we employ an adia¬ 
batic index 7 = 4/3, a ratio of shock to stellar radius with¬ 
out heating ro/r* = 0.4, and an upstream Mach number 
Mo = 5. Examples of the type of initia l con ditio n obtaine d 
with this model can be found in iFernandez fc Thompsonl 
d2009alb Throughout the paper, we make use of a general 
unit system that reflects the dimensionless character of the 
problem. These units and their characteristic values in the 
core-collapse supernova problem are described in Table [Tj 


2.2 Numerical Setup 

We use FLASH3.2 dFrvxell et al.ll2000l : iDubev et al.ll2009l ~) 
to solve the Euler equations in spherical polar coordinates 
(r, 9, </>), subject to the gravity of a point mass and neutrino 
energy source terms, 


U + (v.V)v 

de in t p dp 

df p 2 df 


0 




Qv 


(4) 

(5) 

( 6 ) 


Table 1. Frequently-used quantities and reference values ap¬ 
plicable to core-collapse supernovae. The use of a gamma-law 
equation of state and an arbitrary normalization for the cool¬ 
ing function allow a dimensionless formulation of the problem. 


Symbol 

Quantity 

CCSN Reference value 

r 0 

Initial shock radius 
with B = 0 

100 rioo km 

B 

Heating amplitude 

Equation j2j 

Mns 

NS mass 

1.3 M1.3 .Mq 

M 

Accretion rate 

0.2Mo.2 Mq s -1 

vs 

(2GM ns /r 0 ) 1 / 2 

10 9 - 77 M*g 2 rjjy 2 cm s -1 

to 

rot vs 

E7r%lM~l /2 ms 

Eq 

MroVff 

10 49 - 37 Mo. 2 r 1 1 /oM 1 1/ 3 2 erg 


where p, v, ei n t are the density, velocity, and internal energy, 
and d /At = d/dt + v • V. The system of equations is closed 
with an ideal gas equation of state, and the point mass M ns 
remains constant in time. 

The public version of FLASH has been modified to en¬ 
able the split Piecewise Parabolic Method (PPM) solver to 
operate in 3D spherical coordinates. A detailed description 
of the changes to the code and verification tests of the im¬ 
plementation are presented in Appendix s The code has 
previously been modi fied to en able use of a non-uniform 
grid, as described in [Fernande j (i2012h . 

The computational domain covers the full range of polar 
angles, 9 £ [0, 7 r] and <j) £ [0, 27r], In the radial direction, 
the domain extends from r = r* = 0.4ro to r = 7ro . The 
radial grid is logarithmically spaced with 640 cells, resulting 
in a constant fractional spacing A r/r = 0.45% or 0.26°. 
The meridional grid has constant spacing in cos#, and the 
azimuthal grid is uniform. This angular spacing yields cells 
that subtend a constant solid angle at fixed radius. 

We adopt two resolutions for the angular grid in or¬ 
der to test the sensitivity of the results to this variable. 
The baseline resolution employs 56 cells in the 9 direc¬ 
tion and 192 cells in ij>, for an effective resolution A 9 ~ 
A (j> ~ 2° at the equator. In addition, we evolve two mod¬ 
els at twice the resolution in both 9 and <j> (112 and 384 
grid points, respectively). Our radial resolution is among 
the highest in published 3D hydrodynamic supernova stud¬ 
ies, wh ile our baseli n e ang u lar resolution is compara ble to 
that of lHanke et alJ (120131) . iMezzacappa et akl (l2015fl . and 
iMelson et al.l I 2015blb Our highest angular resolution (1° at 
the equa tor) is only a factor of 2 lower than the highest reso - 
d2012|), ICouch fc O’Connoil (l2014h . 
(2014). 


lution of [Burrows et al 


and 1 Abdikamalov et al. 


The boundary condition in azimuth is periodic. In the 
radial direction, we use a reflecting inner boundary at r = r*, 
and set the ghost cells to the steady-state solution at the 
outer boundary. In the 9 direction, we fill the ghost cells 
with information from active cells located across the axis 
(Appendix 0; for comparison, we also run a model with a 
reflecting boundary condition in 9. 
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4 Rodrigo Fernandez 


2.3 Models Evolved 

We evolve two sequences of models, each with increasing 
heating rate B (eq. El), as shown in Table [2] The two se¬ 
quences differ in the magnitude of the nuclear dissociation 
e at the shock, which changes the strength of the shock 
jump_and t hus the magnitude of the postshock velocity (e.g., 
iThompsonlboOCi ). This difference places each sequence well 
above and below the bifurcation in parameter space that 
determines whether the SASI or convection dominate the 
dynamics. 

The relative importance of each instability is set by the 
relation between the advection rate and convective growth 
rate in the initial condit ion. This can be qua ntified with the 
‘convection’ parameter fFoglizzo et al.ll2006f ) 

x-f TTW (7) 

Jrg pH 

where cubv is the Brunt-Vaisala frequency, and the integral 
runs from the gain radius r g to the shock radius r s . The tran¬ 
sition from SASI- to convection-dominance occurs around 
X — 3. Our SASI-dominated sequence (models names start¬ 
ing with ‘s’) has e = 0, yielding values of x ^ 3 in all initial 
models, while the convection-dominated sequence (model 
names starting with ‘c’) has e = 0.3GM ns /ro, yielding \ 3 
for all values of B. 

For both sequences, we evolve models in ID, 2D, and 
3D around the critical value of B for which an explosion is 
obtained. For ID and 2D, these critical values were obtained 
in Paper I. Here we recompute the same models to maintain 
consistency with 3D in the numerical implementation and 
grid resolutioif]. The critical heating rates for ID and 2D 
models are consistent with the previous value. 

Two types of initial perturbation are applied. In most 
cases, random velocity perturbations in all coordinate di¬ 
rections are applied over the entire computational domain 
at t = 0, with magnitude 0.1% of the radial velocity. In 
a few models, an overdense spherical shell with an i = 
1 angular dependence i s add ed t o the flow (as in, e.g., 
Blond in fc Mezzacappal 120061 or iFernandez fc Thompsonl 
2009bl b The dipole axis of this perturbation is not aligned 
with any cartesian coordinate direction, pointing instead to¬ 
wards 0 = (j) = 45°. This perturbation excites a clean slosh¬ 
ing SASI mode, thus imposing a deterministic initial condi¬ 
tion for probing the effects of dimensionality and resolution. 
Appendix E shows how this type of perturbation is used as 
a test of the ability of the code to maintain different sym¬ 
metries. 

A model is considered to have exploded when the shock 
radius reaches or approaches the outer radial computational 
boundary. For the ID case with finite nuclear dissociation, 
we assign an explosion time at the moment when the peri¬ 
odic oscillation cycle is broken. Explosion times are listed in 
Table [2] Most models are evolved up to a time t = 500to or 
until they explode, whichever comes firsfl 

1 The models of Paper I did not employ a hybrid Riemann solver 
at the shock, whereas our 3D implementation requires such a 
hybrid approach to minimize numerical noise (Appendix [All. 

2 The high-resolution model S08Ll-hr is interrupted earlier than 

500to when it appears that it will not explode, to save computing 
time. 


Table 2. Model parameters and explosion time. Columns 
show from left to right: model name, nuclear dissociation en¬ 
ergy at the shock, heating normalization (eqns. [J-[ 2 ]), initial 
X parameter (eq. [7]), type of perturbation (5v: random ve¬ 
locity, 1=1: overdense shell), and explosion time. A model is 
considered to have exploded when the shock hits or approaches 
the outer simulation boundary. For ID cases with e ^ 0, the 
model ‘explodes’ when the shock oscillation is broken by tran¬ 
sient expansion (e.g. model cld-08). 


Model 

£ 

(GM ns /r 0 ) 

B 

(io- 1 2 3 ) 

X 

Pert. 

texp 

(io) 

ID: 






sld-09 

0 

9 

1.3 



sld -10 


10 

1.6 



sld -11 


11 

1.8 


314 

cld-07 

0.3 

7 

9.1 



cld-08 


8 

11 


350 

2D: 






s 2 d -00 

0 

0 

0 

(5v 


s2d-06 


6 

0.6 



s2d-08 


8 

1.1 

i = i 


s2d-09 


9 

1.3 

(5v 


s2d-09hr 






s 2 d -10 


10 

1.6 


238 

s 2 d- 10 hr 





256 

c2d-06 

0.3 

6 

7.1 



c2d-06hr 






c2d-065 


6.5 

8.0 



c2d-065hr 






c2d-07 


7 

9.1 


423 

c2d-07hr 





339 

3D: 






SOOdv 

0 

0 

0 

5v 


S06dv 


6 

0.6 



S07dv 


7 

0.8 



S08dv 


8 

1.1 


306 

S09dv 


9 

1.3 


204 

SlOdv 


10 

1.6 


219 

S08L1 

0 

8 

1.1 

i = i 

172 

S08Ll-ref 





270 

S08Ll-hr 






S09dv-hr 


9 

1.3 

5v 

187 

C06dv 

0.3 

6 

7.1 

(5v 


C065dv 


7 

8.0 



C07dv 


7 

9.1 


336 

C07dv-hr 





313 


3 RESULTS 

3.1 Transition to Explosion: Overview 

To illustrate the general behavior of SASI-dominated sys¬ 
tems, we describe the evolution of the marginally-exploding 
3D model S08dv (Table [2j) . Initially, a dominant sloshing 
SASI mode grows out of random initial velocity perturba¬ 
tions, as shown in Figure [T}x. As this mode achieves satura¬ 
tion, sloshing modes in orthogonal directions also grow and 
become non-linear, slightly out-of-phase, but the resulting 
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Figure 1. Snapshots in the evolution of the marginally exploding SASI-dominated model (S08dv, top) and convection-dominated model 
(C07dv, bottom), both with random initial velocity perturbations. Each panel shows a slice of the entropy (eq. 0) on the xy plane; to 
ease comparison, the entropy in model C07dv is normalized to the initial postshock value of model S08dv. The SASI-dominated model 
initially develops a sloshing mode (panel a) which saturates (panel b). The subsequent development of a spiral mode (panel c) generates 
a large-scale, high-entropy bubble which leads to the final runaway (panel d). The convection-dominated model initially develops small 
convective plumes (panel e) which gradually grow into larger structures (panels f-g), until a single dominant bubble remains (panel h). 
Note the different spatial scale in panels (d), (g), and (h). 


S08dv (SASI-dominated) 


C07dv (Convection-dominated) 



Figure 2. Explosion geometry near the end of the simulation for the marginally exploding SASI-dominated model (S08dv, left) and 
convection-dominated model (C07dv, right), corresponding to panels (d) and (h) of Figure^ A representative entropy isosurface is shown 
in red, the shock is the outer light grey surface, and the neutron star is the dark grey sphere (the scale is different in the two panels). 
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spiral modc[3 survives only for a fraction of an oscillation 
period and the shock expansion dies down. 

After this lull in shock expansion around t ~ 200to, a 
coherent spiral mode arises, growing quickly in magnitude 
for three full oscillations. By this time, a large high-entropy 
bubble has developed, and the triple-point intersects this 
bubble, interrupting the oscillation. Runaway does not fol¬ 
low immediately, however, with the bubble lingering for a 
period ~ 30to, which allows turbulence to partially break 
the bubble down into smaller parts. This accounts for the 
multiple plumes seen in Figure [TJl, which become more ap¬ 
parent when visualizing the final explosion geometry in Fig¬ 
ure^ If a realistic EOS that includes the energy release from 
alpha particle recombination had been used, the explosion 
pattern would have likely frozen at a smaller radius, and the 
geometry would consist of a single dominant bubble. 

The development of SASI activity is an optimal way 
to rapidly seed large-scale high-entropy bubbles. Whenever 
the shock expands relative to its equilibrium position, the 
Mach number in the frame of the shock increases, and the 
entropy jump correspondingly also increases (e.g., Paper I). 
The for mation of su ch large-scale bubbles in 3D was also 
seen bv lHanke et al.l (I2013T ) when the SASI was active. 

The interplay between bubbles and SASI oscillations is 
also present in 2D (Paper I). In that case, it is found that 
large bubbles cut off accretion to the cooling region, which 
in turn breaks the SASI oscillation cycle. If bubbles are able 
to survive for a time longer than approximately a few SASI 
oscillations, a runaway expansion of the shock will ensue. 
In 3D, the same dynamics appears to take place, but now 
bubbles can have larger sizes given the lack of an axisym- 
metry constraint. As in in 2D, failure to achieve runaway 
expansion in 3D appears to be related to the disruption of 
bubbles due to either accretion plumes and/or turbulence. 
We elaborate on these processes in 1 )3.21 and 1 )3.31 

The dynamics of convection-dominated explosions is 
very different. At early times, numerous small convective 
plumes grow out of the random initial perturbations and 
fill the gain region, as shown in Figure [T| The added con¬ 
vective stresses cause the shock to gradually expand, with 
convective plumes growing in size. This process continues 
until there is a single, dominant convective bubble that ap¬ 
proaches the edge of the computational domain. 

The increase in the shock radius due to convective 
stres s es has been studied extensively (e.g., iMurphv et al.l 
120131: ICouch fe Ottl l2015l l. The formation of a dominant 
structure as the limiting case of bubble growth and/or con¬ 
solidation is an int rinsic proper ty of the non-linear Rayleigh- 
Taylor instability dShartjl 19841 '). A characteristic bubble size 
that subtend s a so lid angle_~ 1 radian was predicted the¬ 
oretically by iThompsonl (I2000h as a result of a balance of 
buoyancy and ra m pressure on t he bu bble, and has bee n 
further studied bv lDolence et all ll2013ll an d|C ouc hi (l2013ll . 

The slow growth of this bubble is in part an artifact of 
the simplified nuclear dissociation prescription, which causes 


3 I Fernand ez! (|201(J) showed that spiral modes can be described as 
a superposition of linearly independent sloshing modes (quantified 
by the real spherical harmonic coefficients of the shock) which are 
out of phase relative to each other. This phase difference does not 
need to be tt/2. 


an excess energy loss relative to t he gravitational potential 
energy as the shock moves out (e.g., F ernandez & ; T homp son 
l2009ah . With a more complete EOS, runaway would have 
ensued at a smaller radius and hence the convective pat¬ 
tern would have frozen before reaching this stage ©• Nev¬ 
ertheless, it is interesting to note that this geometry, with 
a one-sided bubble, has also been obser ved in a simulations 
that include the full equation of state (D olence et al.l 120131 : 
iLentz et al.ll2015h . 

3.2 Dependence on Dimensionality and 
Resolution 

The most interesting result of this investigation is that SASI- 
dominated explosions can take place at lower neutrino lumi¬ 
nosities with increasing dimensionality. Figure [3] illustrates 
the magnitude of the effect, showing the evolution of the 
shock radius with time for systems around the critical heat¬ 
ing rate. Models in 3D can explode at neutrino luminosities 
~ 20% lower than in 2D. 

In contrast, convection-dominated systems show a much 
smaller difference in the critical luminosity between 2D and 
3D within the spacing in B studied here (~ 8%). The un¬ 
derlying difference will manifest itself at smaller intervals of 
B, but at that point the result will also depend on other 
factors such as the numerical resoluti on or the detailed 
form of pre-collapse p erturbations (e.g.. ICouch &; Ottll2015l : 
iMiiller fc Jankall2015h . 

Most of the work that has analyzed the dependence 
of 2D-3D differences in convection-dominated systems has 
concluded that higher resolution is detrimental for 3D (e.g., 
lHanke et al.ll2012 : ICouchll2013l ; I'fakiwaki et al-lfeoMl ). This 
trend has been attributed to the direction in which ki- 
netic energy flows in the turbulent cascades in 2D and 3D 
(lHanke et al.|[2012h . In our models, simply doubling the an¬ 
gular resolution does not prove to be very informative. The 
marginal 2D model at high resolution (c2d-07hr) explodes 
significantly earlier, as expected, but the explosion time of 
the marginal 3D model (C07dv-hr) is not very different- 
The non-exploding 3D model with the highest heating rate 
(C065dv) appears to be starting an explosion on a long 
timescale, but due to the lack of a clear runaway we con¬ 
sider it a failure. Given the significant body of knowledge 
on this type of explosions, we will not concern ourselves fur¬ 
ther with convection-dominated systems. 

The existence of a decrease in the critical luminosity 
with increasing dimensionality for SASI-dominated systems 
is not very sensitive to resolution, as inferred from the fact 
that the model with B = 0.009 explodes with standard 
(S09dv) and high angular resolution (S09dv-hr), at times 
within ~ 10% of each other, and with a heating rate that 
is already below the critical value in 2D. Nevertheless, the 
magnitude of this improvement appears to decrease with 
resolution, as found from the marginally-exploding model 
S08L1. When doubling the angular resolution in 9 and <j>, the 
model does not explode, even though very large shock ex¬ 
cursions still occur (Figure|U). If the full EOS was employed, 

4 Small variations in the explosion time are likely affected by 
stochastic fluctuations, and should not be taken as a primary 
indicator for susceptibility to explosion. 
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Figure 3. Evolution of the shock radius for SASI-dominated models (left) and convection-dominated models (right) with random initial 
velocity perturbations. Top, middle, and bottom panels show models around the critical transition to explosion in ID, 2D, and 3D, 
respectively. For 2D and 3D models, solid lines denote the angle-averaged shock radius, and the shaded area marks the region between 
the minimum and maximum shock radii. Note that 3D SASI-dominated models can explode at a lower neutrino luminosity than in 2D, 
and that non-exploding SASI-dominated models undergo much larger shock excursions in 3D than in 2D ( jpyp . 


this model would likely have exploded due to the additional 
energy from alpha particle recombination at large radius, 
but within the approximations made in this study we con¬ 
sider it a failure. SASI-dominated systems with heating rates 
slightly below the critical value display much larger shock 
excursions in 3D than in 2D, a result that could potentially 
lead to further improvement in explosion conditions once 
the full EOS is included (c.f. H3.3I) . Finally, the decrease in 
the critical luminosity is independent on the boundary con¬ 
dition employed at the polar axis: model S08dv-ref, which 
employs a reflecting boundary condition, explodes at a time 
similar to the fiducial model S08dv. 

We now explore quantitatively the effects of dimension¬ 
ality and resolution in SASI-dominated models. As a base¬ 
line of comparison, we take the marginally-exploding model 
in 3D for which an initial £ = 1 sloshing SASI mode is ex¬ 
cited via an overdense shell (S08L1). This deterministic ini¬ 
tial perturbation allows a straightforward comparison with 
an identical model in 2D (s2d-08Ll), and another in 3D with 
twice the resolution in both 6 and <j> (S08Ll-hr). 

A number of time-dependent diagnostics are employed 
in the analysis. Aside from the angle-averaged shock radius, 


we compute the real l = 1 spherical harmonic coefficients in 
three orthogonal (Cartesian) directions: 

Oi(t) **JdnYi{0,4>)r.(9,<l>,t), (8) 

where r s is the instantaneous shock surface, obtained as a 
pressure-gradient-weighted average for smoothness, and the 
real spherical harmonics basis functions are 


/ 3 

Y{ a ',, y ,z} — a / — {sin 0 cos 0, sin 6 sin </>, cos0} . (9) 

The resulting coefficients can be combined linearly in order 
to measure sloshing modes along arbitrary axeqj. Sloshing 


5 In the 3D models studied here, a sloshing mode is initially ex¬ 
cited along 0 = (p = 45°. Re-defining this direction as the new z 
axis, the real spherical harmonic coefficients transform according 
to: 

a z — y 0.5(ttx + dy) + o>z/y/ 2, 

CLx —y 0.5(tta; T dy') — G-z/\/2, 

CLy —y (CLy — dx) / y/ 2 . 
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Figure 4. Evolution of the shock radius for two SASI-dominated 
models that differ only in their angular resolution. Solid lines show 
the angle-averaged shock radius, with the shaded area marking 
the space between minimum and maximum shock radius. In both 
cases, an (. = 1 sloshing SASI mode is initially excited at an angle 
to the coordinate axes (6 = <f> = 45°). Model S08L1 is just above 
the critical heating rate for explosion. Doubling the resolution 
in both 9 and <f> (model S08Ll-hr) leads to a failure within the 
simulated time. 


SASI modes manifest as sinusoidal oscillations in these co¬ 
efficients. Multiple coefficients oscillating in phase indicate 
sloshing modes along directions not aligned with the coor¬ 
dinate axes. Out-of-phase oscillations between diff erent co¬ 
efficie nts are the signature of spiral SASI modes dFernanded 
l201fll ). The phase difference does not have to be n/2 for the 
composite mode to display the characteristic properties of 
a spiral mode (i.e., angular momentum redistribution). We 
have ignored modes with i ^ 2 given that £ = 1 is expected 
to be the mos t unsta ble osci llat ory SAS I mode at any heat¬ 
ing rate (e.g.. iFhrnandez fc Thompsonll2009bl : Paper I) and 
to keep the analysis concise. 

We diagnose the presence of bubbles with enhanced en¬ 
tropy in the post-shock domain by computing the instanta¬ 
neous fraction of the postshock volume V occupied by fluid 
with entropy higher than a suitably chosen value so (Paper 

I): 


fv 




( 10 ) 


The creation, growth, and destruction of bubbles can be 
characterized quantitatively by following the evolution of 
fv evaluated at multiple values of so- 

In addition, we compute the kinetic energies in the fluc¬ 
tuating components of the flow 


E ± 

E r 



(Vr) 2 ) , 


( 11 ) 

( 12 ) 


where v± is any non-radial velocity component (9 and (j> 
when defined relative to the z axis, or poloidal and toroidal 
otherwise), the mean radial velocity ( v r ) is computed as an 


instantaneous angle-average fe.g.. [Hanke et ahll20ld) . and 
the inner radius for integration ri n is the gain radiuqj. 

Figure [5] shows the time-dependent diagnostics for the 
2D and 3D comparison models (s2d-08Ll and S08L1, respec¬ 
tively). The evolution of the angle-averaged shock radius is 
identical up to time t = 60to, showing the growth and sat¬ 
uration of the sloshing mode that is initially excited. At 
this time, a sloshing mode in an orthogonal direction (a x ) 
has achieved a noticeable amplitude. Thereafter, the aver¬ 
age shock radius undergoes a larger shock excursion in 3D 
relative to 2D. Nonetheless, this initial expansion dies down 
over a time equivalent to several SASI oscillations. At time 
t = 115to a new expansion begins, and the 3D model transi¬ 
tions to a runaway. Meanwhile, the 2D model does not show 
any significant increase in its average shock radius, achieving 
a statistical steady-state. 

The spherical harmonic coefficients indicate that the z 
component, initially excited, achieves the same saturation 
amplitude in 2D and 3D. This is e xpecte d from estim ates 
of the effect of parasitic instabilities (IGuilet et alj|2010ll and 
also manifests in the SASI-only tests in Appendix [X] The 
final phase of spiral mode growth that precedes explosion in 
3D is characterized by a comparable amplitude in all three 
(out-of-phase) real spherical harmonic coefficients. Interest¬ 
ingly, in this second expansion phase none of the orthogonal 
sloshing components in 3D achieves a higher saturation am¬ 
plitude than the 2D mode except after runaway expansion 
has set in. 

Also shown in Figure [5] are the non-radial kinetic ener¬ 
gies in the gain region. Initially, the poloidal kinetic energy 
in 3D (defined relative to the sloshing axis excited) is iden¬ 
tical to the 9 kinetic energy in 2D. The energies begin to 
evolve differently once the a x and a y modes are excited, 
manifesting as an increase in the toroidal kinetic energy. 
This toroidal energy decreases as the poloidal energy sat¬ 
urates, at the time when the shock is retracting. The sum 
of poloidal and toroidal energies remains comparable to the 
9 energy in 2D. Once the spiral mode achieves noticeable 
amplitude at t ~ 100to, both toroidal and poloidal energies 
increase monotonically. The time at which the final expan¬ 
sion of the shock begins coincides with the time when the 
sum of the non-radial kinetic energies in 3D exceeds that in 
2D. By the end of the exploding simulation, this difference 
is a factor of several. 

The volume fraction diagnostic shows that in fact the 
2D model achieves larger overall fractions of the volume with 
high entropy than the 3D model. This would normally indi¬ 
cate that the 2D model achieves more favorable conditions 
for explosion. Upon closer inspection, however, we find that 
while the 2D model contains structures with higher entropy 
than in the 3D model, these structures occupy disconnected 
volumes, as shown by the xz entropy slices in Figure [5Ji-b. 
These structures also live for a shorter time, as indicated by 
the oscillatory pattern in the fv diagnostic for the 2D model. 
Thus, achieving explosion conditions is not only related to 
the size of the bubbles or the magnitude of the entropy in 


6 When comparing models with and without heating ( 33.3D . the 
inner radius for integration is taken to be that where the angle- 
averaged sound speed is maximal; c.f. Paper I. 
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Figure 6. Entropy slices in SASI-dominated models with the 
same parameters and initial conditions except dimensionality 
and/or resolution, as labeled. Panels (a-b): Slices on the xz plane 
for 3D and 2D models. In the former, the plane is defined rela¬ 
tive to the sloshing mode that is initially excited (0 = 4> — 45°). 
While the 2D model reaches higher entropy, the sizes of the bub¬ 
bles are smaller than in 3D. Panels(c-d): Slices on the xy plane 
in 3D models with different resolution. The plane is perpendicu¬ 
lar to the sloshing axis. The model with high angular resolution 
displays a larger degree of small-scale turbulence. Compare with 
Figs. \E\ and [7] 


Figure 5. Comparison between two SASI-dominated models that 
differ only in their dimensionality: S08L1 (3D, marginally explod¬ 
ing) and s2d-08Ll (2D, non-exploding). In both cases, an £ = 1 
sloshing mode is initially excited; for the 3D model, the sloshing 
axis is along 6 = (p = 45°. Panel (a): Average shock radius. Panel 
(b): Real spherical harmonic coefficients (eq. [H]) defined relative 
to the initial sloshing axis in 3D or along the symmetry axis in 
2D. Panel (c): Non-radial kinetic energies. In 3D, the poloidal and 
toroidal directions are defined relative to the initial sloshing axis. 
Panel (d): Fraction of the post-shock volume that has entropy 
higher than a given value so (eq. 0), as labeled. The solid and 
dashed lines correspond to the 3D and 2D, models, respectively. 


these structures, but also involves the time during which 
these structures can survive. 

Regarding the effect of resolution, Figure [7] shows the 
diagnostic quantities for the two 3D models that differ only 
in their angular resolution, one exploding (S08L1) and the 
other non-exploding (S08Ll-hr; c.f. Figure[4]). The evolution 
of both models is qualitatively very similar until t ~ 130fo, 
after which the evolution bifurcates. The spherical harmonic 
coefficients show that while the evolution of the sloshing 
mode that was initially excited is very similar in both cases, 


the amplitude of a x is larger in the lower resolution model at 
the time when explosion begins. In other words, a stronger 
spiral SASI mode can be the difference between explosion 
and failure. 

This subtle difference is also reflected in the kinetic en¬ 
ergies in the gain region. Both radial and non-radial ener¬ 
gies are initially larger in the higher resolution model, as ex¬ 
pected from the known de pendence of the SASI grow th rates 
with angular resolution ((Fernandez fe Thompsonl l2009bh . 
The energies in the two models remain close to each other 
during the period of shock retraction, and bifurcate only 
around the time when runaway sets in. 

The volume fraction diagnostic fv is more informative 
in this case, indicating that the lower resolution model is 
able to maintain a larger volume with high-entropy relative 
to the high-resolution model after time t ~ 70fo- 

Figure [S] also shows two snapshots of the entropy on 
the xy plane (defined relative to the initially excited slosh¬ 
ing mode) at similar stages prior to the bifurcation in the 
evolution of these two models. The snapshots suggest that 
the smaller scale turbulence in the higher resolution model is 
related to more efficient fragmentation of high-entropy bub- 
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Figure 7. Comparison between two SASI-dominated models that 
differ only in their angular resolution: S08L1 (standard resolu¬ 
tion, marginally exploding) and S08Ll-hr (high resolution, non- 
explodig). In both cases, a sloshing mode along 0 = <f> = 45° is 
initially excited. Panel (a): Average shock radius. Panels (b) and 
(c): Real spherical harmonic coefficients (eq.[8| defined relative to 
the initial sloshing axis. Panel(d): Transverse and radial kinetic 
energies (eq. 1111 and 1121) . Panel (e): Fraction of the post-shock 
volume that has entropy higher than a given value so (eq. W), 
as labeled, with the entropy is defined as in equation {3J). Solid 
and dashed lines correspond to standard and high resolution, re¬ 
spectively. 


bles. We interpret the failure to explode in model S08Ll-hr 
as due to the shorter life of these high-entropy bubbles given 
the more efficient disruption by turbulence. This interpreta¬ 
tion is a conjecture at present, requiring further study to be 
established as the definitive cause for the failure of higher 
resolution models. An interesting question is whether this 
shredding of bubbles can be characteriz ed in terms of para - 
sitic instabilities as for SASI modes (e.g., [Guilet et alll201 Oh . 



Figure 8. Relation between shock expansion and kinetic energy 
in selected non-exploding SASI-dominated models, illustrating 
the effect of neutrino heating and dimensionality (models s2d- 
00, s2d-06, SOOdv, S06dv; see Table [2j. Panel (a): Average shock 
radius. Panel (b) and (c): Non-radial kinetic energy (eq. 1111 ) for 
models with B=0 and B=0.006, respectively. Panel (d): Radial 
kinetic energy in the fluctuating component of the flow (eq. |12| ). 
The operation of spiral modes in 3D leads to a larger amount of 
non-radial kinetic energy than in 2D. Episodic shock excursions 
are related to the formation of large-scale high-entropy bubbles 
when neutrino heating is included. When the bubbles are shred¬ 
ded and an explosion fails to start, the shock retreats and the 
cycle resets. 


3.3 Spiral mode with heating: large shock 
excursions 

Another interesting result of this investigation is that SASI- 
dominated models in 3D that do not explode, but which 
are close to the critical neutrino luminosity, display shock 
excursions that are much larger than in their 2D counter- 
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parts. This effect is visible in Figure [3] for the model with 
B = 0.007. Here we investigate the origin of this effect, and 
discuss the potential for a larger difference between 2D and 
3D once a realistic EOS, with alpha particle recombination, 
is included in the evolution. 

For this analysis, we carry out two additional non¬ 
exploding models, each in 2D and 3D: one with no heat¬ 
ing (s2d-00 and SOOdv, respectively) and another with B = 
0.006 (s2d-06 and S06dv, respectively). Figure [8^i compares 
the evolution of the average shock radius for these four mod¬ 
els. Two features stand out. First, the 2D models have sys¬ 
tematically smaller shock radii than their 3D counterparts. 
Second, models with heating undergo episodes in which the 
shock executes large excursions. In the 2D model with heat¬ 
ing, these excursions last for a shorter time and have smaller 
amplitudes than in 3D. 

Inspection of the 3D model with heating shows that 
these excursions are associated with the appearance of large- 
scale, high-entropy bubbles. An expansion of the shock 
against a supersonic incident flow leads to an increase in 
the Mach number in the frame of the shock, leading to a 
larger entropy jump and hence to the formation of a high- 
entropy bubble with a size comparable to the portion of the 
shock that is expanding (e.g., Paper I). As is the case in 
2D, the formation of large enough bubble leads to the in¬ 
terruption of the SASI cycle. If a runaway expansion fails 
to ensue, the bubble is shredded, either by turbulence or by 
external bulk motions (e.g., downflows), and the shock re¬ 
tracts. A larger shock radius means a larger volume subject 
to neutrino heating where liigh-entropy bubbles can grow. 

Figure [8jo-c also shows the non-radial kinetic energies 
for the four comparison systems (c.f. eq. in]). Given that 
models without heating are included in the comparison, the 
lower radius for computation of these energies is that where 
the angle-averaged sound speed is maximal (e.g., Paper I). 
The resulting (total) non-radial kinetic energies are larger 
in 3D than in 2D, whether neutrino heating is included or 
not. It is not surprising then that the average shock radius 
is larger in 3D than in 2D, as more kinetic energy inside 
the shock provides additional support against t h e ram pres¬ 
sure of the collapsi ng stellar core (e.g.. iMuruhv et aLll2013l : 

Icouch fc otdEoish . 

For models without heating, in which non-radial mo¬ 
tions are solely driven by the SASI, the non-radial kinetic en¬ 
ergies in the fully non-linear phase satisfy Eg D ~ E# ~ Eg D . 
Hence, the spiral modes of the SASI are able to approx¬ 
imately double the total non-radial kinetic energy in the 
post-shock domain. Inspection of Figure [5jo as well as tests 
in Appendix |A] indicate that the saturation amplitude of in¬ 
dividual sloshing modes, quantified by the real spherical har¬ 
monic coefficients (eq. [9]), is very similar in 2D and 3D_[this 
is also a prediction from the analysis of lGuilet et al.l[201 Cf). 
Recall that two or three linearly-independent real spheri¬ 
cal harmonic coefficients oscillating in phase are indicative 
of a single sloshing mode, thus their individual amplitudes 
will not reach the maximum possible value (their root-mean- 
square sum will). A spiral mode appears to be able to tap 
into the full amplitude of its individual constituent (orthog¬ 
onal) sloshing modes, hence roughly doubling the kinetic 
energy. 

For the cases with heating, the non-radial kinetic energy 
in 2D is slightly larger than the individual 8 or <j> energies 


in 3D, but still smaller than their sum. The episodic shock 
excursions in 3D are concurrent with increases in the non- 
radial kinetic energy. However, the overall magnitude of the 
non-radial 3D energies are very close to those of the model 
without heating. 

The key feature that sets apart the model with large 
shock excursions from others is the magnitude of the fluctu¬ 
ating radial kinetic energy (eq. 1121 1. as shown in Figure [8jl. 
Whereas in 2D the radial kinetic energy is comparable to 
that in the 8 direction whether heating is included or not, 
in 3D the total non-radial kinetic energy exceeds its radial 
counterpart by about ~ 50% during most of the time. How¬ 
ever, during episodes of transient shock expansion, the ra¬ 
dial kinetic energy of the model with B = 0.006 shows large 
increases that match or even exceed the total non-radial ki¬ 
netic energy. 

We interpret this increase in radial kinetic energy as 
being a consequence of the large high-entropy bubbles hav¬ 
ing internal structure and being convective, as can be seen 
from Figure [l] In a fluid where buoy ant convecti on oper¬ 
ates, the relation E r ~ Eg + E$ holds dMurphv et al.ll2013l ). 
The fact that the radial kinetic energy remains large only 
for a short period of time is related to the short lifetime of 
the large-scale bubbles, which are shredded prior to shock 
retraction. 

The large shock excursions shown by the 3D SASI- 
dominated models with neutrino luminosities close to criti¬ 
cal hold promise for further reductions in this critical heat¬ 
ing rate once a more realistic EOS is included. If the ef¬ 
fect of alpha particle recombination is included, the ra¬ 
tio of nuclear dissociation energy to gravitational poten¬ 
tial energy at th e s hock d oes not increase with radius (e.g. 
IFernandez fe Thompsonll2009al l. and hence the system loses 
less energy upon expansion. A characteristic radius r a where 
this effect becomes dynamically important is that were the 
nuclear binding energy of alpha particles equals the gravita¬ 
tional binding energy, 

r a — 254Mi. 3 km. (13) 

Taking rg = 100 km for a characteristic stalled shock radius 
without neutrino heating, shock excursions such as those of 
model SOOdv should very likely lead to a runaway expansion. 
This conjecture will be tested in future work. 


4 SUMMARY AND DISCUSSION 

We have investigated the properties of SASI- and 
convection-dominated core-collapse supernova explosions 
in three dimensions using parameterized hydrodynamic 
simulations. This approach allows to isolate the effect of 
hydrodynamics instabilities on the explosion mechanism 
from uncertainties in progenitor models, the EOS of dense 
matter, or an incomplete treatment of neutrino effects. By 
generating sequences of models well within the parameter 
regimes where either the SASI or convection dominates, we 
characterized the effect of adding a 3rd spatial dimension 
on the dynamics. Our main results are the following: 

1. - SASI-dominated explosions do exist in 3D, and 
their evolution is qualitatively different from convection- 
dominated explosions. (Fig. [I]). For systems in which both 
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the mass accretion rate and protoneutron star radius change 
slowly relative to the thermal time in the gain region, the 
value of the convection parameter x (eq. [I]) at the time 
when the shock stalls is a good predictor of the explosion 
path to be taken by the system, if the neutrino luminosity 
is high enough. 

2. - SASI-dominated systems can explode with a lower 
neutrino luminosity in 3D than in 2D (Fig. [3], Table [5]). 
This difference is related to the ability of spiral modes 
to generate more non-radial kinetic energy than a single 
sloshing mode, resulting in a larger average shock radius 
and hence generating more favorable conditions for the for¬ 
mation of large-scale, high-entropy bubbles (Fig.0. For the 
baseline angular resolution employed here (AS ~ Ac/ ~ 2° 
at the equator), the critical neutrino luminosity in 3D is 
~ 20% lower than in 2D. 


3. - Convection-dominated explosions show a much smaller 
change in the critical heating rate with increasing di¬ 
mensionality. In fact, no significant difference was found 
between 2D and 3D for the spacing in heating rate explored 
here (~ 8%). This result is consistent with previous studies 
(e.g.. lHanke et ahll2012 .1. Doubling the angular resolution of 
the marginally exploding models yields earlier explosions in 
2D and minor changes in explosion time for 3D (Table 0. 


4. Doubling the angular resolution in 9 and </ decreases 
the difference in critical neutrino luminosity between 
2D and 3D for SASI-dominated explosions (Figs. [4] and 
0. We interpret this reduction as a consequence of the 
higher efficiency of 3D turbulence at higher resolution for 
disrupting the large-scale, high-entropy bubbles needed to 
launch an explosion (Fig. [6]). The exact dependence of this 
increase in the critical neutrino luminosity with resolution 
in 3D was not obtained, but based on the healthy explosion 
obtained in a 3D model with 10% lower heating rate than 
the marginal 2D model (B = 0.009, both at high resolution; 
Table 0, we infer that this difference is shallower than 
BcSt x (A(9Ac/) -0 ' 085 . This sensitivity to resolution also 
highlights the need for a better understanding of the 
systematic uncertainties associated with the use of finite 
volume methods to model supernova flows, in which the 
spatial resol ution is limited by current computational re - 
sources (e.g., IPorter fe Woodward! 19941 ; [Radice et al.ll2015t l. 


5. - Non-exploding SASI-dominated models have larger 
average shock radii in 3D than in 2D, which we interpret 
to be a consequence of the higher non-radial kinetic energy 
generated by spiral modes (Fig. 0. While the saturation 
of individual sloshing modes appears to be very similar 
in 2D and 3D (Fig. [5j Appendix [At. a spiral mode can 
generate more kinetic energy than a single sloshing mode 
(e.g., Figure IA2I) . This result is in agreement with the 
findings of lHanke et al.l (l2013h . and in conflict with the 
interpretation of Iwakami et al.l ([2008) that saturation in 
3D is smaller because the same amount of kinetic energy 
is shared between a larger number of modes. The exact 
limiting factor to the energy available for a spiral mode was 
not investigated here, but it is worth pursuing in future 
studies. 


6 . - Very large shock excursions can result in 3D when the 
heating rate is close to (but lower than) the critical value for 
an explosion in SASI-dominated models. These excursions 
are simultaneous with the appearance of high-entropy 
bubbles generated by the expanding shock, and are smaller 
in magnitude (by a factor ~ 2) in 2D. The non-radial kinetic 
energy is generally larger than the radial kinetic energy 
when spiral modes operate, except when large bubbles are 
formed and the shock expands, in which case the relation 
E r ~ Eg + E$ holds. We conjecture that this is due to the 
internal structure of bubbles, which is likely convective. 
The large shock excursions open the possibility of further 
improvements in explosion conditions in 3D when the effect 
of alpha particle recombination is included. 


W hile the behavior of the 27 Mq model of lHanke et ali 
(l2013h during the time when spiral modes operate appears 
to be consistent with our results, the period of spiral-SASI- 
driven shock expansion is too short to make an unambiguous 
parallel. The maximum amplitudes of the sloshing modes 
involved and the total transverse kinetic energies exceeded 
the values obtained in 2D for a brief period. Unfortunately, 
just as the shock begins to turn around and grow at a rate 
faster than in 2D, the Si/SiO composition interface is ac¬ 
creted through the shock, triggering a large expansion that 
makes the model convective again. 

The relation between growth of the non-radial ki- 
netic energy and appro ach to explosion was no ted by 
I Murphy fe Meakinl (1201 il l and lHanke et al.l (l2012lf . Intu¬ 
itively, one expects accreted material to spend more time 
in the gain region when the flow is non-radial and the 
gain region is large r due to the larger shock radius (e.g., 
iMurphv fe Burrowsl l2008l i. An interesting way to think 
about the reduction in critical neutrino luminosity when 
going from ID to 2D/3D is that due to the non-laminar 
stresses, a given shock radius requires less thermal pressure 
to be maintained and thus less neutrino heating to drive 
an e xpansion from this position than if the flow were lami¬ 
nar (ICouch fc O tt 2015). This is equivalent to saying that, 
in analogy with accretion disks, the effective gravity felt by 
the fluid inside the shock is smaller due to the action of bulk 
or turbulent stresses, 


5eff 


GM(r) < ct 


(14) 


where n| uct = (v% — ( v r ) 2 + Vq + v^). This is also consis¬ 
tent wit h the cr itical r.m.s. Mach number threshold pro¬ 
posed by I Muller fc Jankal d2015ll . Regardless of the concep¬ 
tual framework used to interpret it, the improved ability of 
spiral modes to extract non-radial kinetic energy from the 
accretion flow appears to improve conditions for explosion. 

While we have found a region of parameter space in 
which 3D is more favorable to explosion than 2D, it is fair 
to ask whether this parame ter space is actuall y realized in 
Nature. The 27 Mq model of lHanke et al.l d2013h is genuinely 
SASI-dominated until the Si/SiO interface is accreted, as di¬ 
agnosed by the evolution of x, which lingers around 2 before 
the composition interface falls throug h the shock. In con¬ 
trast, the same progenitor evolved by ICouch fe O ’Connor! 
[201J) (using /heat = 1.00) yields x closer to 3 when 
the shock stalls, showing that the treatment of neutrino 
transport can be influential in determining the evolution- 
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ary path taken by the system. More generally, our knowl¬ 
edge of two key external ingredients to the supernova mech¬ 
anism is by no means complete. First, progenitor stars are 
likely to be members of int eracting binary systems (e.g., 
ISana et~al . 2012; Smith|2014|) , intrinsically multi dimensional 
Ve.g.. lArnett fc Meakinll201ll : ICouch et al.ll201fll l and with a 
degree of rotatio n yet to be convincingly established (e.g., 
IFulIer et aij 120151 1. Second, the EOS of matter at supra¬ 
nuclear densit ies remains unce rtain despite significant recent 
progress (e.g.. lLattimei)l2012h . 

In the end it may well tu r n out to be that, as sug¬ 
gested by lAbdikamalov et al.l ( 2011 1. the SASI is rele¬ 
vant mostly for high accretion rate progenitors that fail 
to explode. Because neutron star formatio n generall y pre¬ 
cedes collapse to a black hole (e.g.. lO’Connor fc Qttll2011I f. 
the same predictions for the neutrino and gravitational 
wave emission for systems with strong S ASI activity (e.g., 
iMiiller et al.l 120131 : iTamborra et al.l I2014T ) should apply. In 
fact, the higher accretion rates in failed systems are ex¬ 
pected to yield a more i ntense n eutr ino signal than ex¬ 
ploding models (lO’Connor fc Ot.TI 120131 ) . This should im¬ 
prove prospects for identific ation of a Galactic ev ent even 
if dedicated surveys (e.g., iKochanek et al.l 120081 ') fail to 
detect t he electromagne t ic signature predicted for t hese 
systems iNadezhinj 1980l; Lovegrove fc Wooslevl 120131 ; IPirol 


1 20131 : iKashivama fe Quataertl 2015ll . 


If on the other hand progenitor properties turn out 
to favor SASI-dominated explosions once uncertainties have 
been resolved, the effect of pre-collapse perturbations on the 
development of a spiral SASI mode remains as a potential 
obstacle. The advective-acoustic cycle involves the interplay 
of coherent perturbations that require a relatively smooth 
background flow to develop. The work of lGuilet et alj (120101 1 
indicates that large amplitude perturbations can disrupt the 
acoustic feedback in the cy cle and lead to in e ffectiv e growth. 
Furthermore, the work of IMiiller fc Jankal (l2015h suggests 
that if physically-motivated pre-collapse perturbations are 
used, the shock acquires a nearly static asphericity, with 
detrimental conditions for the development of the SASI (but 
providing an alternative channel for the generation of non- 
radial kinetic energy). 

Our results prove the principle that SASI-dominated 
explosions can provide better conditions for explosion in 3D 
than in 2D. For the effect to be considered robust, however, 
a key set of physical ingredients needs to be added. Future 
work will address the response of these types of explosion 
to the inclusion of a realistic EOS (including nuclear recom¬ 
bination), rotation in the accretion flow, the feedback from 
the accretion luminosity on the gain region, and the time- 
dependence of the accretion rate and neutron star radius. 
These studies will help clarify issues such as the behavior of 
the 20 M p) model with sophistic ated neutrino transport re¬ 
ported in lTamborra et all (|2014l l. which in 3D undergoes an 
extended phase of SASI activity but fails to explode, while 
succeeding in 2D (T. Janka, private communication). 

A more realistic (but still parametric) study would 
also allow exploration of issues such as the elementary 
processes that mediate the transition to explosion (e.g., 


Pe jcha fc Thompsonll20l3 : lFernande3l2012l : IMiiller fe Jankal 
20151 1. the amount of angular momentum imparted to the 


neutron star by spiral modes, w hich depends on the mass 
cut at the time of explosion (e.g.. lGuilet fc Fernandea l2014f 


as well as t he potential for spir al modes to amplify magnetic 
fields fe.g.. lEndeve et ahll2012l l in an explosion context. The 
results of such studies would inform the analysis of more 
sophisticated models, in which the high computational cost 
precludes a systematic exploration of parameter space. 

Ultimately, the magnitude of the improvement in effi¬ 
ciency on the neutrino mechanism introduced by the SASI in 
3D, if all conditions are favorable, is of the order of ~ 10%. 
While this contribution is relatively modest, it can help to 
tilt systems that would be otherwise marginal towards a ro¬ 
bust explosion (as is the case when, e.g., the neutrino cross 
secti ons are modified at t he ~ 10% level in the right direc¬ 
tion; iMelson et al .1 [20 1 Sail . 
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APPENDIX A: EXTENSION OF THE SPLIT 
PPM SOLVER IN FLASH3.2 TO 3D SPHERICAL 
COORDINATES 

Here we describe the extension to 3D spherical polar coor¬ 
dinates of the split PPM solver present in the public version 
of FLASH3.2. We then describe tests of this modified code. 
The test cases can be used to diagnose any hydrodynamics 
code in 3D spherical coordinates. 


Al Implementation 


The split PPM so lver em bedded in FLASH3.2 is based on 
PROMET HEUS (iFrvxell et alll989h, which implements the 
method of lColella fc Woodward! " 19841) . The public FLASH 
version requires only a few modifications to run on a grid 
defined in 3D spherical coordinates. The modifications are: 


(i) Allowing 3D spherical geometry in the subroutine 
Hydro_detectShock, and computing the appropriate veloc¬ 
ity divergence; 

(ii) Modifying the ratio of time step to cell spacing dt/Ax 
in the subroutine hydro_ld to account for the proper cell 
spacing in the azimuthal direction, dx$ = r sin 6d(j> 

(iii) Implementing a safeguard in the subroutine avisco, 
whenever division by sin# leads to singular behavior in the 
velocity divergence. This is the only subroutine where this 
divisor is evaluated at the cell face, where it can vanish if 
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the grid extends to the polar axis. In practice, a small floor 
value of sin 8 (~ 10 -6 ) can be imposed. 

(iv) Modifying the arguments to the subroutine 
Hydro_detectShock, called by hy_ppm_sweep in order 
to use the correct coordinates. 

In addition to these basic modifications, we have exper¬ 
imented with two additional changes that can improve the 
accuracy of the code: 

(i) The fictitious accelerations that account for the cen¬ 
trifugal and Coriolis forces in curvilinear coordinates are 
computed with velocities evaluated at the beginning of the 
time-step. For consistency with PPM, these accelerations 
should be time-centered. We implement this modification 
by evaluating 1/2 of the acceleration forward in time after 
the PPM update. This is straightforward because the force 
along a given direction depends linearly on the correspond¬ 
ing velocity component. Tests indicate that this modification 
leads to minor differences in the results. 

(ii) The default reflecting boundary condition in the 8 
direction at the axis can be improved by making use of the 
data across the axis to fill the ghost cells. For example, for 
the first active cell next to the axis (with center coordinate 
8 = 9i), the contents of its neighboring ghost cell can be set 
to satisfy 

A(r, —0\,4>) = A(r, 9i,<f> + n), (Al) 

where A is an arbitrary variable. The sign of the transverse 
velocities ve and must be changed to account for the sign 
change in the corresponding unit vectors across the axis. 
This type of operation requires ghost cell exchange between 
processors that contain cells located at opposite sides of the 
axis. 

Additional changes are needed only when using a non- 
uniform grid. These changes involve generating the non- 
uniform coordinates, using the exact cell volumes and areas 
(instead of linearized approximations), and passing the cor¬ 
rect vector of cell spacings to the PPM routines instead of 
a constant val ue. T hese modifications are further described 
in lFernande3 (120121 s ) . 

A2 Tests 

We test the reliability of the hydrodynamic solver in 3D 
by using the parameterized accretion shock setup described 
in H2.2I without neutrino heating. The system is then only 
unstable to the SASI, with no convection. By exciting indi¬ 
vidual modes of system, the numerical solution can be com¬ 
pared directly with predi ctions f rom linear stability analysis 
((Fernandez fc Thompsonll2009bl l. 

The parameters of the setup are an initial ratio of star 
to shock radius r*/ro = 0.5, adiabatic index 7 = 4/3, no nu¬ 
clear dissociation at the shock (e = 0 ), and the same param¬ 
eterized neutrino source term in equation © but without 
neutrino he ating (B = 0). This is the same set of parame¬ 
ters used in lFernandej (|201ol ). which has the advantage that 
only the fundamental i = 1 SASI mode is unstable. In the 
absence of rotation, the modes are degenerate in m, serv¬ 
ing as a good diagnostic of the isotropy of the code. This 
is a global problem involving subsonic flow, where a deli¬ 
cate balance between pressure gradients and gravity needs 


Table Al. Test models evolved. Columns from left to right 
show model name, type of initial perturbation applied, use of 
a hybrid Riemann solver at shocks, type of boundary condi¬ 
tion used at the polar axis in the 8 direction, and root-mean- 
square amplitude of the l = 1 spherical harmonic coefficient 
along the excited dipole axis (for reference, the 2D value is 
Aai/ro = 0.544). The initial perturbation has the form of an 
overdense shell with angular dependence given by a real spher¬ 
ical harmonic: Y z , Y x , or Y,i = Y z /y/2 + (Y x + Y y )/2. When no 
explicit perturbation is applied, initial transients result in a 
spherically-symmetric perturbation. For comparison, ID and 
2D versions of each model are also evolved. 


Model 

Pert. 

Hybrid? 

Axis Bnd. 

Aai/r 0 

T-L0-hyb 

none 

yes 

reflect 


T-LO-std 


no 



T-Llz-ref 

Y z 

yes 

reflect 

0.534 

T-Llx-ref 

Y x 



0.530 

T-Lld-ref 

Y d 



0.544 

T-Llz-trm 

Y z 

yes 

transmit 

0.533 

T-Llx-trm 

Y x 



0.526 

T-Lld-trm 

Yd 



0.552 


to be maintained. Any obvious errors in the code manifest 
immediately. 

The tests differ in the type of initial perturbations ap¬ 
plied and in the boundary condition at the axis. SASI modes 
are excited by dropping an overdense shell with a given an¬ 
gular dependence, as described in H2.21 For comparison, ID 
and 2D versions are also evolved. We used the standard an¬ 
gular resolution (A 9 ~ A <j> ~ 2° at the equator, H2.2I) for 
all tests. Table [All summarizes the parameters of each test 
simulation. Model names start with T for “test”, and then 
indicate the type of initial perturbation and boundary condi¬ 
tion employed (e.g., T-Llz-ref corresponds to 1= l sloshing 
mode along the z axis, with reflecting boundary condition 
at the axis). 

The tests probe four aspects of the code: (1) the ability 
to maintain spherical symmetry and agreement with ID and 
2D versions, (2) the ability to maintain axisymmetry and 
its agreement with the 2D version, (3) the effect of using a 
reflecting or ‘transmitting’ boundary condition at the axis, 
and (4) the isotropy of the code in 3D (symmetry along an 
arbitrary axis). 

A 2.1 Spherical symmetry 

The first test simply lets the accretion flow evolve in the 
absence of perturbations. Initial transients generate an (. = 0 
perturbation that damps on a timescale of several 100 to, 
because the corresponding SASI mode is stable. 

To achieve a clean test setup, it is important to damp 
numerical perturbations that arise spontaneously. We find 
that significant numerical noise arises from the shock in 3D 
models if the standard PPM solver is used in all cells. This 
is illustrated in Figure [All which shows the angle-averaged 
shock radius for models T-LO-std and T-LO-hyb, as well as 
the magnitude of the non-radial kinetic energies in the 9 and 
(j> directions. 

It is well known that a hybrid Riemann solver can elim¬ 
inate numerical problems at shocks aligned with the grid, 
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Figure Al. Test of the ability of the code to maintain spherical 
symmetry in the absence of explicit perturbations. Panel (a): Av¬ 
erage shock radius as a function of time for ID, 2D, and 3D models 
without initial perturbations. Also shown is a 3D model that does 
not use the hybrid Riemann solver. Panel (b): Transverse kinetic 
energies for 3D models without perturbations, showing hybrid 
Riemann (black) and standard split PPM (red). Using the hybrid 
Riemann solver keeps spurious non-radial velocity perturbations 
from growing too much above the numerical noise. 

particularly the odd-even dec oupling instability and carbun¬ 
cle phenomenon (|Quir kj[l994]). FLASH3.2 allows the use of a 
HLLE solver inside shocks to increase dissipation and damp 
these numerical instabilities. Figure lATl shows that when this 
hybrid Riemann solver is used, spurious numerical perturba¬ 
tions (quantified by the non-radial kinetic energies) remain 
close to the numerical noise, growing slowly. 

Using the hybrid Riemann solver also enables excel¬ 
lent agreement between ID, 2D, and 3D. This agreement 
is reached only above a certain radial resolution, however, 
due to the steep density gradient at the base of the flow, 
which results in the irregular misidentification of shocks at 
low resolution. The baseline radial resolution adopted in the 
study (a fractional radial spacing Ar/r ~ 0.45% or 0.26°) 
is the lowest resolution for which this agreement between 
different dimensions is reached. 


A 2.2 Symmetry around the z axis 

The next test probes the ability of the code to maintain 
axisymmetry around the 2 axis. We evolve a 3D model for 
which an i = 1 sloshing mode is excited along the z axis 



Figure A2. Test of the ability of the code to maintain axisym¬ 
metry. Panel (a): Real spherical harmonic coefficients (eq. [8]) 
for model T-Llz-ref and an identical run in 2D. In both cases an 
1=1 SASI sloshing mode along the z axis is initially excited. 
The code remains axisymmetric until the SASI saturates. Panel 
(b): Kinetic energies in the polar (6) and azimuthal (</>) directions 
for model T-Llz-ref. Upon saturation of the SASI, the azimuthal 
kinetic energy grows rapidly from numerical noise and triggers 
the growth of modes orthogonal to the z axis. 


(T-Llz-ref). Given our choice of parameters, only the funda¬ 
mental l = 1 mode should be unstable and therefore a clean 
sinusoidal oscillation of the corresponding l = 1 spherical 
harmonic coefficient (eq. [8] ) should be obtained. 

Figure fA~2k shows the evolution of the 2 real spherical 
harmonic coefficient (eq. [H]) for model T-Llz-ref, and com¬ 
pares it with a 2D version that excites an t = 1 mode in the 
same way. The 3D model remains axisymmetric up to the 
moment when the sloshing SASI mode saturates at t ~ 60fo- 
Thereafter, the azimuthal kinetic energy experiences rapid 
growth out of numerical noise, as shown in Figure IX2b . lead¬ 
ing to the excitation of transverse sloshing modes in the x 
and y directions. This growth in transverse motion origi¬ 
nates in numerical noise at the polar axis. Nonetheless, the 
sloshing mode along the 2 direction remains nearly identical 
between 2D and 3D up to a time t ~ 70fo, or ~ 24, 000 time 
steps. At later times the primary difference lies in the phase 
of the oscillation. 

The rms fluctuation in the a z coefficient between t = 75 
and 150to is Aai/ro = 0.544 and 0.534 in 2D and 3D, respec¬ 
tively, with a relative difference of 2%) in favor of 2D. Also, 
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time [t 0 ] 


Figure A3. Non-radial kinetic energies - defined relative to the 
ic-axis (eq. lA2l and lA3t - as a function of time for models T-Llx- 
ref (black lines) and T-Llx-trm (red lines), in which a sloshing 
SASI model along the m-axis is initially excited, and which use 
reflecting and transmitting boundary conditions in the 9 direction 
at the axis, respectively. Solid lines show poloidal and dotted lines 
toroidal kinetic energy. 


note that the difference in phase between the real spherical 
harmonic coefficients of the 3D model shows that a spiral 
mode is triggered. 


A2.3 Axis boundary condition 

To obtain a quantitative measure of the improvement gained 
by using a ‘transmitting’ boundary condition in 9 (eq. IA1I) 
instead of a reflecting z axis, we evolve two models in which 
an i = 1 sloshing mode along the x axis is excited. 

As a diagnostic, we define poloidal and toroidal veloci¬ 
ties relative to the rr-axis: 


xy x z g 

- Vy H- V Z - V X 

(A2) 

rg rg r 

z y 

Vy + Vz 5 

(A3) 


e Q 


with q = sjy 1 + z 2 . These velocities are the analog of vg and 
v^, respectively. In theory, the sloshing mode should remain 
symmetric around the x axis, with the toroidal velocity Vtor 
serving as a diagnostic for numerical errors. 

Figure IA3I shows the evolution of the poloidal and 
toroidal kinetic energies in models T-Llx-ref and T-Llx- 
trm, constructed with the corresponding velocities defined 
in equation and (IA3I) . In contrast to the sloshing mode 
along the z axis (Fig. IA2I) . for which the azimuthal kinetic 
energy undergoes rapid growth only after the primary SASI 
mode has saturated, the toroidal kinetic energy in both mod¬ 
els with an x sloshing SASI mode grows exponentially from 
the beginning. This illustrates the larger amount of numer¬ 
ical noise experience by this mode given the grid geometry. 

Using the transmitting boundary condition results in 
a lower amplitude of the toroidal kinetic energy at early 
times, by a factor of several, relative to using a reflecting 
boundary condition. Nonetheless, both types of boundary 
condition lead to very similar results in the evolution of the 
primary mode. The lower level of noise motivates the use 



time [t Q ] 


Figure A4. Test of the isotropy of the code. Curves show the real 
spherical harmonic coefficient (eq. [8]) corresponding to the axis 
along which a sloshing SASI mode is initially excited: z (solid 
gray), x (dashed red), and a diagonal axis (z/V2 + [x + y]/2 , 
dotted blue). The top panel shows models that use a reflecting 
boundary condition in the 6 direction at the polar axis, and the 
bottom panel shows models that employ a transmitting boundary 
condition ( LV2.3l i. 

of the transmitting boundary condition as the default for 
production runs. 

A2-4 Symmetry along an arbitrary axis: isotropy 

The final test diagnoses the isotropy of the code, which we 
assess by comparing the evolution of sloshing SASI modes 
excited along different directions: (1) the z axis, which is 
the natural symmetry axis of the grid given the coordinate 
system, (2) the x axis, which while orthogonal to the z axis 
is still aligned with a coordinate line on the equatorial plane, 
and (3) a diagonal axis defined by 

d=-kz+i(* + y). (A4) 

This direction lies along 6 = <j> = 45°, and hence it is not 
aligned with any cartesian coordinate direction. 

Figure [ATI shows the real 1=1 spherical harmonic co¬ 
efficients (eq. [8] ) along the corresponding sloshing axis for 
the last six models of Table [ATI The evolution is identical 
in all cases up to the point where the primary SASI mode 
saturates, demonstrating the correctness of the implemen¬ 
tation. 
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Models that use a reflecting boundary condition at 
the axis show mutual agreement in the non-linear phase 
for a longer time relative to models with a ‘transmitting’ 
axis. This agreement extends particularly to the maximum 
amplitude of the sloshing mode. In contrast, models with 
‘transmitting’ boundary condition show closer agreement in 
the oscillation phase at late times, but the amplitudes are 
not identical. The root-mean-square fluctuation of the real 
spherical harmonic coefficient, computed between t = 75 
and 150£o is shown in Table [All All values are within a few 
percent of each other, and within a few percent of the 2D 
value (0.544). 
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